Dynamical Transition in Interaction Quenches of the One-Dimensional Hubbard 

Model 
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We show that the non-equilibrium time-evolution after interaction quenches in the one dimen- 
sional, integrable Hubbard model exhibits a dynamical transition in the half-filled case. This transi- 
tion ceases to exist upon doping. Our study is based on systematically extended equations of motion. 
Thus it is controlled for small and moderate times; no relaxation effects are neglected. Remarkable 
similarities to the quench dynamics in the infinite dimensional Hubbard model are found suggesting 
dynamical transitions to be a general feature of quenches in such models. 
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I. INTRODUCTION 

In the last years, seminal experimental setups have 
been developed combining a very good decoupling of 
the quantum systems from their environment with a 
high degree of controllability of the system's parame- 
ters. This renders the observation of the temporal evo- 
lution of closed quantum systems for long times possi- 
ble. Switching the internal parameters provides tools 
to investigate systems out of equilibrium. In optical 
lattices, the internal parameters such as hopping and 
particle-particle interaction can be manipulated^—. Fur- 
thermore, various pumb-probe experiments based on ul- 
trafast spectroscopy have been developed^. 

These developments have triggered extensive theoret- 
ical studies of physics far from equilibrium, based on a 
large variety of analytical and numerical tools^— . The 
goal is to qualitatively understand and to quantitatively 
describe the dynamics of quantum systems far from ther- 
mal equilibrium. How do such states evolve? Do they 
relax towards equilibrium? How does this happen and 
on which time scales? 

These issues are relevant in all dimensions. But so 
far the infinite dimensional and the one-dimensional 
(ID) case have attracted the greatest interest. Mainly, 
this is due to the possibility of approximation-free re- 
sults in both cases. The infinite dimensional case 
(d —> oo) is amenable by dynamical mean-field the- 
ory (DMFT)i£~— for which it is exact, also in the non- 
equilibrium case^. The ID case is amenable to analyti- 
cal (e.g. bosonizatio n 7 ' 13 ! 19 : 20 ) and numerical (e.g. time- 
dependent density-matrix renormalizatio n 9 ' 10 ! 21 : 22 ) ap- 
proaches. In addition, there is a fascinating conceptual 
issue inherent to ID systems: All extended non-trivial 
integrable quantum systems are ID and the existence of 
an exhaustive number of conserved quantities influences 
the dynamics strongly because it implies a macroscopic 
number of conserved quantitities^. Still, relaxation in 
local quantities due to dephasing may occur—. 

One efficient way to realize states far from equilib- 
rium are interaction quenches^. The interaction value 



is changed abruptly increasing it 7 .i 9 . ~ 11 i 13 i 19 or decreas- 
ing . We focus on interactions which are suddenly 
switched on with Fermi seas as initial states. This sim- 
plifies the evaluation because all correlations are known 
analytically, but it is not essential for the applicability 
of our approach which is based on extended equations of 
motion^ 3 -. 

For an understanding of the dynamics of quenched sys- 
tems also the behavior on short and moderate times di- 
rectly after the quench matters. Relaxation, if it hap- 
pens, may set in directly or coherent oscillations may 
take place before relaxation occurs on a longer time scale. 
The occurrence of intermediate time scales, on which the 
system does not yet reflect the long-time behavior, is sub- 
sumed under the term "prethermalization" n i 25 ] 26 . These 
issues are particularly difficult to study in strongly cor- 
related systems. 

Interestingly, evaluating the DMFT equations after in- 
teraction quenches of the half-filled system two regimes 
occurred^ which are characterized by qualitative differ- 
ences in the time dependence of the jump An(t) in the 
momentum distribution at the Fermi level. For small val- 
ues of the interaction U, this quantity decreases gradually 
from its initial value of unity corresponding to the non- 
interacting Fermi sea and a prethermalization plateau 
can be discerned. For large values of U, the jump An(t) 
is characterized by strong oscillations. The key differ- 
ence is that An(t) displays zeros in the strong-interaction 
regime while they are absent for weak interaction. 

This observation of two qualitatively different dynam- 
ics could be retrieved in a variational Gutzwiller ap- 
proach, evaluated in Gutzwiller approximation, by Schiro 
and Fabrizio^ 7 .. Their treatment maps the quantum dy- 
namics to a classical mechanics problem. True relaxation, 
however, is neglected in this way. But the resulting equa- 
tions allow for an analytic evaluation. The oscillatory 
regime for weak quenches and the one for strong quenches 
are separated by a singularity indicating a dynamical 
transition. Note that the Gutzwiller approximation be- 
comes exact for infinite coordination numbe r 28 : 29 so that 
the similarity between the DMFT and the Gutzwiller re- 
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suits may not suprisc. 

In view of these results, it is our goal to show that a 
Hubbard model with completely different properties dis- 
plays essentially the same dynamical transition. Thus 
we study the one-dimensional Hubbard model. It is 
different from the Hubbard model examined by DMFT 
and Gutzwiller techniques in two important aspects: 
(i) Scattering between the excitations is controlled by 
momentum conservation in 1D22t— while momentum 
conservation is suppressed at internal vertices in large 
dimension s 1 6 ' 1 7 . (ii) The macroscopic number of con- 
served quantities in the integrable ID Hubbard model^ 
should constrain the non-equilibrium dynamics so that 
it differs significantly from what is seen in higher dimen- 
sions. Yet we find that the dynamical transition is also 
present in ID with even quantitative similarities. 

The article is set up as follows. Next, the model, 
the quench, and our approach to them are introduced. 
In Sect. Ill, the results for the time-dependent jump in 
the momentum distribution are presented. In Sect. IV, 
we analyse the data with respect to the occurrence of a 
dynamical transition. Finally, conclusions are drawn in 
Sect. V. 



II. QUENCHED HUBBARD MODEL 

For the quench, we start from a Fermi sea and switch 
on the interaction abruptly. Thus the Hubbard Hamilto- 
nian becomes time-dependent and reads 

# Hu = - j ( 6 Ls> + h - c -) + u <t) (!) 

with a local repulsion U(t), where (dj) create (anni- 
hilate) a particle with spin a at site j and hj = CjCj. 
We study U(t) = Q(t)U > and define the band width 
W = 4 J as the natural energy scale. 

The approach used is a systematically controlled ex- 
pansion of the Heisenberg equations of motion for cj CT — . 
By commuting the interacting H after the quench re- 
cursively with Cj CT , i.e., by applying the Liouvillian, we 
obtain differential equations for the prefactors of the ex- 
pansion of £j a (t) in more and more monomials of c\ a and 

Cj at t = 0. We perform this calculation in real space. 
The application of one commutation is a 'loop'. Roughly, 
each loop multiplies the number of tracked monomials by 
a factor 3. We are able to realize up to 11 loops. In Fig.Q] 
the convergence of time-dependent results with the loop 
number is shown. Our approach bears similarities to re- 
cent calculations based on recursively constructed Hilbert 
spaces^. But we stress that our approach is operator- 
oriented. The evaluation of expectation values is done 
only at the end at the time instant of interest. 

The evaluation of the differential equations on the one- 
loop level describes the time dependence of the initial 
operator cj a exactly in linear order in t. By each loop 



we increment the depth of the hierarchy of the equations 
of motions by one and thus the time dependence of the 
initial operator cj a is precisely captured up to order i" 
for n loops. Finally, we solve the differential equations 
numerically so that also higher powers of t are gener- 
ated. The conceptual asset of the approach is that it 
works directly on the infinite lattice by exploiting trans- 
lational invariance. No relaxation effects occurring up to 
the considered order are neglected. A quantitative analy- 
sis of the convergence of the approach and of its accuracy 
is presented in App. [5] 



III. TIME-DEPENDENT JUMP IN THE 
MOMENTUM DISTRIBUTION 




t [1/W] 

FIG. 1: (color online) Upper panel: Jump An(t) for the 
half-filled Hubbard model for various loop numbers. Lower 
panel: An(t) for increasing U (from top to bottom at small t: 
U/W = 0.125, 0.25, 0.5, 0.8, 1.0, 1.25, 1.5, 2.0, 5.0) in 11 loops. 



We focus on the momentum distribution rik(t) := 
(c' k a c k a )(t) where k is the wave vector—. In particular, 
we study the jump An(t) := ?ifc F +o(0 — n-k F -o{t) at the 
Fermi vector k-p with Art(0) = 1 for the initial Fermi sea. 
After the quench, An shows slow relaxation to zero and 
damped oscillations. The upper panel of Fig. [T] depicts 
the jumps for increasing number of loops for the half- 
filled Hubbard model. Good convergence is obtained for 
11 loops up to about t w 10/ W . The precise value up 
to which the data is reliable depends on the details, see 
App. [5] The lower panel of Fig. [1] displays the data for 
various values of U . 

The slow relaxation seen in Fig. [T] is character- 
istic for a ID model as can be understood from 
results for Tomonaga-Luttinger models obtained by 
bosonizatio n 7 ' 13 i 19 ' 20 . It is found that power laws instead 
of exponential relaxation occur. Our data agrees with 
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this expectation see data in Ref. |3J 
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FIG. 2: (color online) Jump An(t) for the half-filled Hubbard 
model in ID (black, solid lines) and on the Bethe lattice with 
infinite coordination number (gray, dashed lines from Ref.—). 
Dotted line: Second order result An = 1 — U 2 t 2 /A. 



Our first key observation stems from the comparison of 
the ID data with exact DMFT data for a Bethe lattice^. 
The common lore expects a crucially distinct behavior 
due to the differing dimensionality and the integrability 
of the ID Hubbard model^. The scattering in ID is 
strongly restricted due to momentum conservation^ - — 
while this conservation is irrelevant at internal vertices 
in DMFT-i&i£. Yet Fig. [2] shows qualitatively similar re- 
sults for larger values of U. For U ;$ W quantitative dif- 
ferences prevail beyond t = 2/W so that it is difficult to 
discern qualitative similarities. But for interaction val- 
ues beyond the band width W coherent oscillations in 
An(t) occur— which agree in that the minima appear at 
almost the same values. To stress that this is not simply 
an effect of the leading order in t the second order result 
An(t) = 1 - U 2 n(2 - n)t 2 /A (n filling factor) is included 
in Fig. [2] The amplitudes of the oscillations, however, 
differ significantly. We conclude that the periods of the 
oscillations are a candidate for similarities while the am- 
plitudes and their decay are more strongly depending on 
the model. In ID, a slow power law relaxation is generic 
while exponential relaxation is expected in higher dimen- 
sions in agreement with the DMFT data. These striking 
similarities indicate that for large U the physics is gov- 
erned by local processes as is explained below. 



IV. DYNAMICAL TRANSITION 

We return to the oscillations oc cos(27ri/T) which ap- 
pear both for strong quenches and for weak quenches. 
Schiro and Fabrizio derived an analytic formula for the 
period T of the oscillations at half-filling 



where K(x) is the complete elliptic integral of the first 
kind and U c the value where the Mott transition occurs in 
the Gutzwiller approach^!. It is given by SE^n where 
-Ekin is the kinetic energy of the half-filled Fermi sea. 
Away from half-filling, the logarithmic singularity in ([2]) 
is smeared out. 

What does our data yield for the oscillation periods? 
In absence of relaxation, thus without changes in the os- 
cillatory amplitudes, it is straightforward to read off the 
period T. In our data, this is less obvious. To determine 
the period T we fit a straight line to the data such that 
it is a double tangent to two points close to the first two 
minima, see inset of Fig. |4j From the first instant we 
deduce the period T = 2to. The double tangent is con- 
structed to take into account that the oscillatory part 
may be small and sit on top of a non-oscillatory relax- 
ation. In case only one minimum can be determined, we 
read off its instant in time to find to. 

Representative results are depicted in Fig. @] for vari- 
ous fillings. For comparison, the analytic formula © is 
plotted as well. At and around half- filling (n = 1) and 
for interaction values U smaller than about U c /2 the os- 
cillation period T hardly depends on U. The rather flat 
curve in this regime is in agreement with the semiclas- 
sical finding for half-filling @ although the value of T 
differs roughly by a factor of two which may not surprise 
in view of the very different models. Note that T stays 
finite for U — > 0. Its value in this limit depends on the 
filling. Data for other fillings can be found in App. [Bland 
in Figs. [3] and 2] below. We stress that all oscillations in 
the weak quench regime do no show zeros in An(t) so 
that An(t) > 0. 

For quenches to large U and in particular for U — > oo 
the curves for all fillings converge well to the asymptotic 
behavior T = 2ir/U. This value corresponds to Rabi os- 
cillations in the local two-level system which is given by 
a local singly and doubly occupied site. In other words, 
the lattice behaves as if it were made up from indepen- 
dent Hubbard atoms in first approximation. This finding 
stems from the basic two-loop calculation leading to 

4,f(*) = h (-l,t) : cL lit : +h Q (0,t) : c^-j- : 

+ VI, t) : 4 it : +hi(0, 0, 0, t) : 4,f4,J,^)4 : ( 3 ) 

where site can be any site on the chain because of trans- 
lational invariancc. The coefficients follow the differential 
equations 



Ti / Ti\ 

U- (0,0,0,i) 



T 



^K{2U/U C ) for 

4° 
V,. 



±K{U C /2U) for 



U <U C 
U > U c , 



(2) 



d t h o {0,t) = -Jih (-l,t) - Jih (l,t) 

(4a) 

d t h (-l,t) = -Jiho(0,t) (4b) 
d t h (l,t) = -Jih o {0,t) (4c) 
d t hi(0, 0, 0, t) = Uih o (0, t) + (1 - n)hi(0, 0, 0, t) (4d) 

for the coefficients. This result covers the leading order in 
t independent of the value of U because the commutation 
of the interaction part in (JTJ) with the local operators, i.c, 



4 



those at site 0, in ((3j) does not yield any other operator 
term than those that we included. In this sense, the ap- 
proach used becomes exact for U —> oo. The surrounding 
lattice sites act as damping bath. 

At half-filling, the above differential equations can be 
solved analytically and we find for the oscillation period 
T = 2n/ v / U 2 + W 2 /2. For finite doping the numeri- 
cal solution is easily done. In both cases, T = 2ir/U is 
the leading order result in an expansion in 1/U. The 
equations of motion approach is well-controlled for large 
values of U because the number of operators generated 
by commutation with the interaction part does not grow 
infinitely. Their number is finite for a fixed number of 
sites involved because the corresponding Hilbert space is 
finite. The subleading order W/U 2 is exactly captured 
if all operators acting on two sites are included which 
happens in the seven-loop calculation. Clearly, local pro- 
cesses dominate for large values of U . 

Returning to the two regimes of quenches, we stress 
that only the oscillations at half-filling display zeros in 
An(t) after the strong quenches. For any finite doping 
(n / 1) the oscillations persist, but do not reach zero 
any more in their minima as shown in detail in Fig. [3] 
Even 2% of doping arc sufficient to shift the first mini- 
mum from zero upwards to a finite value within numerical 
accuracy, see right panel in Fig. [3] 



of a singular behavior which may be a jump (indicated 
by a dashed line) or a logarithmic singularity (the peri- 
ods rise by about a factor two above their U — > value) . 
The similarity to the semiclassical result is surprising and 
indicates that indeed the two regimes of quenches are sep- 
arated by a dynamical transition. The small shift of the 
position of the anomaly from U = U c /2 to U fa 0A3U C is 
of minor importance. 
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FIG. 3: (color online) Oscillations in An(t) for var- 
ious dopings. In the left panel, the dopings are 
from top to bottom at the first minimum: n = 
0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0. In the right panel, 
the dopings are from top to bottom at the minimum: n = 
0.8, 0.82, 0.84, 0.86, 0.88, 0.9, 0.92, 0.94, 0.96, 0.98, 1.0. Clearly, 
zeros occur only at half-filling. 

The regimes of weak and of strong quenches are sep- 
arated by an anomaly catching the eye in Fig. U at 
U ~ U c /2. From our data we cannot identify a singu- 
larity because there are some uncertainties in the deter- 
mination of T, see App. [C] for error bars and further 
discussion, and we cannot make statement about times 
beyond ss 15/W. Nevertheless, our finding is suggestive 



FIG. 4: (color online) Period of the oscillations in An(t) as 
function of U (U c = SW/n in ID); n denotes the total filling 
factor per site. The dotted curve depicts the semiclassical 
Gutzwiller result dSpl. 



We emphasize that all oscillations after strong 
quenches, to the right of the dashed part of the half- 
filled curve, display zeros An(t) while no zeros are found 
for weak quenches, to the left of the dashed part. This 
holds, however, only at half-filling. For finite doping, 
even strong quenches do not imply zeros so that the tran- 
sition ceases to exist for any doping, see Fig. [3] Thus the 
data constitutes very strong evidence for a qualitative, 
dynamical transition between both regimes at half-filling. 
At finite doping, it appears to be washed out immedi- 
ately. The anomaly, however, disappears gradually: We 
see in Fig. [4] that about 20% of doping are required to 
make the anomaly vanish completely, cf. App. [B] 

The finding agrees with the semiclassical, relaxation- 
free Gutzwiller result in infinite dimensions^ and with 
the available DMFT datai^. Since such similar behavior 
occurs in extreme cases such as a one-dimensional in- 
tegrable model with important momentum conservation 
and an infinite dimensional model with suppressed mo- 
mentum conservation the conclusion suggests itself that 
the dynamical transition is a general feature in quenched, 
half-filled Hubbard models. 



5 



V. CONCLUSIONS 

Summarizing, we studied interaction quenches in the 
one-dimensional Hubbard model by an approach based 
on equations of motion. It is systematically controlled in 
the depth of the hierarchy; here up to 11 commutations 
were carried out. The accuracy of the data is system- 
atically controlled by the comparison of the results for 
different loop numbers. 

We analyse the timc-dcpcndencc of the jump in the 
momentum distribution starting from the initial Fermi 
sea. Slowly decaying oscillations are found. At half- 
filling, two qualitatively different regimes appear: For 
strong quenches, the jumps display zeros, for weak 
quenches they do not. In the dependence of the oscil- 
lation period T on the interaction U a clear anomaly ap- 
pears separating both regimes. It indicates a dynamical 
transition. Since this scenario is in surprising agreement 
with previous findings in infinite-dimensional model a 12 ! 27 
we suggest that this scenario holds for half-filled Hubb- 
bard models in general. 
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Appendix A: Technicalities 

We are interested in momentum distributions after in- 
teraction quenches. They can be computed by Fourier 
transformation of the one-particle equal time propaga- 
tors 



G(r,t) = (FS|c(r,£) c t (0,i)|FS) 



(Al) 



where the expectation value is taken with respect to the 
noninteracting Fermi sea |FS) which represents the initial 
state before the quench. The r stands for a site on the 
lattice under study, here a chain in one dimension (ID). 
The time-dependent operators c and are represented 
by the following ansatz 



r (r,t) =fl+ (F F ft) +.. 



/ r 



(A2) 



where (L') denote various superpositions of particle 
(hole) creation operators. For instance, the single parti- 
cle creation is given by 



f ?= E Em<u)^ 

|5|gw max t a 



(A3) 



with time-dependent prefactors ho(S,t). The shifts <5 
can be bounded from above by v max t where v max is 
the maximal velocity in the sense of the Lieb-Robinson 



theorei 



,38,39 



Physically, this means that there is a max- 



imum velocity with which the essential effects of the 
Hamiltonian after the quench can travel. Of course, ex- 
ponentially small effects will be neglected. Hence, in 
order to describe the dynamics correctly up a certain 
time one has to include processes up to a certain spatial 
ranged. This idea is also used in other approaches^. 

To calculate the time dependence of the prefactors we 
use the Heiscnbcrg equation dtA(r, t) = i H, A(r, t) for 

the time derivative of an operator A. On calculating the 



commutator 



we encounter two cases. The 

commutation of the noninteracting part of the Hamil- 
tonian Hq leads to a shift of the fermionic operators 
in space, whereas the commutation with the interaction 
term H mt may additionally create or annihilate particle- 
hole pairs ft Lt. Iterating this process then leads to the 
ansatz (|A2|) . With each commutation more terms with 
higher number of particles involved are created. Thus 
the amount of terms grows exponentially: At 11 loops we 
deal with up to 5 • 10 5 monomials in the Hubbard model 
and a set of differential equations with about 2- 10 7 terms 
on the right hand side. 

The differential equations of the prefactors can be 
solved numerically with the initial conditions /io(0, 0) = 1 
and hi(r,t) = Vi ^ 0. Because each commutation com- 
prises one order in time t a calculation with n commu- 
tations provides results for c'(t) which are exact up to 
order t n . We stress, however, that we do not use a series 
expansion in t but solve the (truncated) set of differential 
equations numerically. 

Due to the proliferating number of additional terms 
arising within the calculation we have to restrict ourselves 
to a finite number of commutations. As terms appearing 
during the last commutation for the first time lead to 
an overestimation of the weight loss for the one particle 
terms, we omit them to improve the convergence. A cal- 
culation with m commutations performed in this way is 
called an m-loop calculation. 

To illustrate how each commutation improves the re- 
sult the absolute difference between an m-loop calcula- 
tion and the 11-loop calculation is shown for a quench 
with U = LOW in a double logarithmic plot in Fig. 
[5] Based on this difference a runaway time is defined 
as the time t beyond which the difference takes values 
larger than a certain threshold which we set to 0.01. This 
threshold is depicted in Fig. as dashed line. 

The resulting inverse runaway time is shown in Fig. [6] 
as function of — where m denotes the number of loops 
performed as described above. The curve shows a power 
law increase of the runaway time as function of — , i.e., 
a power law decrease of the inverse runaway time, with 
an exponent of about 1.8. Note that this indicates that 
the convcrcncc is supcrlinear which is favorable for its 
application in practical calculations. 
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FIG. 5: (color online) Absolute difference of the jump An m (t) 
at various numbers of loops m relative to the 11-loop re- 
sult Anii(t) for the half-filled Hubbard model. Dashed line: 
Threshold for the determination of the runaway time. 




0.35 



1/m 



FIG. 6: (color online) Double logarithmic plot of the inverse 
runaway time vs. the inverse number of loops of the corre- 
sponding calculation. Dashed line: Power law fit of the data 
with an exponent of about 1.8. 



Appendix B: Oscillations for Various Filling Factors 

The period of oscillations for various values of the fill- 
ing is shown in Fig. [7] as function of the compactified 
interaction 



(u /2)+u ' interaction U c denotes the val- 
ues of the Hubbard interaction, where the Mott transi- 
tion occurs in the Gutzwiller approximation 2 ^ 29 ' 40 . It 
takes the value 8W/ir in ID. 

Upon doping the anomaly at U ~ U c /2 is gradually 
reduced until it vanishes completely at around n « 0.8, 
leaving behind a monotonic behavior of the period. In 
contrast to the gradual disappearance of the anomaly, 
the zeros for strong quenches to the right side of the 
anomaly vanish immediately upon doping. This was 
illustrated in Fig. [3] 
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FIG. 7: (color online) Period of the oscillations in An(t) as 
function of U (U c = 8W/n in ID); n denotes the total filling 
factor per site. 



Appendix C: Systematic Errors in the 
Determination of the Period 



In the determination of the periods of oscillation sys- 
tematic inaccuracies due to the finite amount of loops 
occur. There arc two main sources: (i) The finite num- 
ber of loops and (ii) the way the period T is extracted 
from the data. The effect of the finite number of loops 
is estimated by comparing results from the calculation 
with 11 loops with calculations with a smaller number of 
loops. This yields estimate (i) for the systematic error. 

In order to estimate the effect of the fitting procedure 
we compare the periods T determined from the double 
tangents described in the main article to an alternative 
analysis. In this alternative procedure we read off the 
positions t' of the first minimum of An(t) and deduce 
V = 2t' . The absolute value \T-T'\ yields estimate (ii) 
for the systematic error. 

Finally, the maximum of both estimates (i) and (ii) is 
used as the error of the period. These errors are shown 
by error bars for fillings n = 1.0, 0.8 and 0.5 for some ex- 
emplary values of the interaction in Fig. [8] In the regime 
of large interactions, the errors are negligible. Around 
the anomaly, they are largest as could be expected. But 
the position and shape of the anomaly arc hardly affected 
by the systematic errors. Thus we are confident that our 
conclusions are valid. 

If a double tangent can be constructed as in the inset 
of Fig. 3 in the main article, a third way to determine 
the period is by T" = t\ — to. This procedure is not 
possible for all points so that we do not use it for the 
determination of the curves T(U). For small values of U, 
the third approach yields periods shorter by up to 1/W. 
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FIG. 8: (color online) Period of the oscillations for the fillings 
n — 1.0,0.8,0.5. The error bars depict the error of the cal- 
culated period for some exemplary values of the interaction. 
They are estimated as described in the main text. 
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